
// Vector/Matrix test routine

// $Id: VectorTest.cpp 2 2012-05-22 17:55:16Z gerry@gaboury.biz $

#include "Vector.hpp"
//#include "VectorStore.hpp"
#include "Matrix.hpp"
#include "MatrixSolver.hpp"
//#include "MatrixSymmetric.hpp"
#include <iostream>
#include <cstdlib>
#include <ctime>

using namespace std;

int
main (int argc, char** argv)
{
  try {

    Vector<double> V1(4,3), V2(4), V3(4);

    cout << "V1 = " << V1 << endl;

    for (int k = 0; k < 4; ++k )
      {
	V1[k] = k;
	V2[k] = 2*k;
      }

    V3 = V1 + V2;
    double prod = V1 * V2;
    cout << "prod = " << prod << endl;

    cout << "V3 = " << V3 << endl;

    Matrix<double> M1(4,4);
    M1.identity();
    cout << "M1 (identity) =\n" << M1;

    if ( M1.isSymmetric() )
      cout << "M1 is symmetric" << endl;
    else
      cout << "M1 is asymmetric" << endl;

    srand48(time(0));

    for ( int i=0; i < 4; ++i )
      {
	M1[i][i] = drand48();
 	for ( int j=0; j < i; ++j )
 	  M1[j][i] = M1[i][j] = drand48();
      }

//     M1[0][0] = 1;
//     M1[1][1] = 2;
//     M1[2][2] = 3;
//     M1[2][0] = 5;
//     M1[3][3] = 4;
//     M1[1][3] = 0.5;

    cout << "M1 = \n" << M1;

    if ( M1.isSymmetric() )
      cout << "M1 is symmetric\n" << endl;
    else
      cout << "M1 is asymmetric\n" << endl;

    //MatrixSymmetric<double> MSym ( 4 );

    Vector<double> V4 = M1 * V3;
  
    cout << "V4 = M1 * V3 = " << V4 << endl;

    Matrix<double> M2 = M1 * 2;

    cout << "M2 = \n" << M2;

    MatrixSolver<double> MS( M1 );

    MS.LU_Decomp();

    cout << "MS (decomposed) =\n" << MS;

    cout << "M1 determinant = " << MS.Det() << endl;

    Vector<double> V5 = MS.LU_Solve ( V2 );

    cout << "V5 = " << V5 << endl;

    MS.Invert();

    cout << "MS (inverted) =\n" << MS;

    if ( MS.isSymmetric() )
      cout << "MS is symmetric\n" << endl;
    else
      cout << "MS is asymmetric\n" << endl;

    Matrix<double> M3 = MS * M2;

    cout << "M3 (should be identity * 2) = \n" << M3;

    Matrix<float> MF1(2,2),MF2(2,2);
    MF1[0][0] = 1;
    MF1[0][1] = 2;
    MF1[1][1] = 4;

    MF2[0][0] = 2;
    MF2[1][0] = 1;
    MF2[1][1] = 3;

    Matrix<float> MF3 = MF1 * MF2;

    MatrixSolver<double> MSym(2,2);
    MSym[0][0] = 2;
    MSym[1][0] = MSym[0][1] = -2;
    MSym[1][1] = 5;

    cout << "before Cholsky decomposition =\n" << MSym << endl;
    MSym.Cholesky_Decomp();
    cout << "after Cholsky decomposition =\n" << MSym << endl;

    Matrix<double> MD1(2,2),MD2(2,2);
    MD1[0][0] = 1;
    MD1[0][1] = 2;
    MD1[1][0] = 3;
    MD1[1][1] = 4;

    MD2[0][0] = 0;
    MD2[0][1] = 5;
    MD2[1][0] = 6;
    MD2[1][1] = 7;

    Matrix<double> KP = MD1.KronProd(MD2);
    cout << "Kronecker Product" << endl;
    cout << KP;

    // generate a symmetic matrix
    Matrix<double> Msym = MD1 * MD1.Tr();
    MatrixSolver<double> Msolve = Msym;
    cout << "symmetrix matrix" << endl;
    cout << Msym.row() << ", " << Msym.col() << endl;
    cout << Msolve << endl;
    Msolve.Cholesky_Decomp();
    Matrix<double> Minv = Msolve.Cholesky_Inverse();
    cout << "inverse" << endl;
    cout << Minv.row() << ", " << Minv.col() << endl;
    cout << Minv << endl;
    cout << "product - should be identity" << endl;
    Matrix<double> Mprod = Msym * Minv;
    cout << Mprod << endl;

    cout << "--> end main routine <--" << endl;
  }
  catch (exception& ex) {
    cerr << ex.what() << endl;
  }
}


